Multistate B- to A- transition in protein-DNA Binding - How well is it described by current AMBER force fields?

When DNA interacts with a protein, its structure often undergoes a significant conformational adaptation, usually involving a transition from B-DNA towards the A-DNA form. This is not a two-state, but rather a multistate transition. The A- and B- forms differ mainly in sugar pucker (north/south) and glycosidic torsion χ (anti/high-anti). The combination of A-like pucker and B-like χ (and vice versa) represents the nature of the intermediate states between the pure A- and B- forms. Here we study how the A/B equilibrium and the A/B intermediate states at protein-DNA interfaces are modeled by current AMBER force fields. Eight diverse protein-DNA complexes and their naked (unbound) DNAs were simulated with OL15 and bsc1 force fields and an experimental combination OL15χOL3. We found that while the geometries of the A-like intermediate states agree well with the native X-ray geometries, their populations (stabilities) are significantly underestimated. Different force fields predict different propensities for A-like states growing in the order OL15 < bsc1 < OL15χOL3, yet all underestimate A-like form populations. Interestingly, the force fields seem to predict the correct sequence-dependent A-form propensity, as they predict larger populations of the A-like form in unbound DNA in those steps that acquire A-like conformations in protein-DNA complexes. The instability of A-like geometries in current force fields significantly alters the geometry of simulated protein-DNA complexes and destabilizes the binding motif, suggesting that refinement is required to improve description of protein-DNA interactions in AMBER force fields.Communicated by Ramaswamy H. Sarma.

Assessing the Current State of Amber Force Field Modifications for DNA─2023 Edition.

Advances in molecular dynamics (MD) software alongside enhanced computational power and hardware have allowed for MD simulations to significantly expand our knowledge of biomolecular structure, dynamics, and interactions. Furthermore, it has allowed for the extension of conformational sampling times from nanoseconds to the microsecond level and beyond. This has not only made convergence of conformational ensembles through comprehensive sampling possible but consequently exposed deficiencies and allowed the community to overcome limitations in the available force fields. The reproducibility and accuracy of the force fields are imperative in order to produce biologically relevant data. The Amber nucleic acid force fields have been used widely since the mid-1980s, and improvement of these force fields has been a community effort with several artifacts revealed, corrected, and reevaluated by various research groups. Here, we focus on the Amber force fields for use with double-stranded DNA and present the assessment of two recently developed force field parameter sets (OL21 and Tumuc1). Extensive MD simulations were performed with six test systems and two different water models. We observe the improvement of OL21 and Tumuc1 compared to previous generations of the Amber DNA force. We did not detect any significant improvement in the performance of Tumuc1 compared to OL21 despite the reparameterization of bonded force field terms in the former; however, we did note discrepancies in Tumuc1 when modeling Z-DNA sequences.

Probabilistic Interpretation of NMR J-Couplings Determines BI-BII State Equilibria in DNA.

Interpretation of 3JP,H3' NMR scalar spin-spin coupling constants in DNA becomes more reliable by including distinct structural states such as BI and BII, using the weighted-static or, better still, the recently implemented adiabatic-MD (Ad-MD) method. The calculation method employs an adiabatic ("Ad") dependence of 3JP,H3' coupling on NMR-assigned torsion angle, ε, weighted by P(ε) probability distribution calculated by molecular dynamics (MD). Ad-MD calculations enable cross-validation of the bsc1, OL15, and OL21 force fields and various parametrizations of the Karplus equation describing the dependence of 3JP,H3' coupling on ε torsion (KE). The mean absolute deviation of Ad-MD 3JP,H3' couplings from the experimental values in Dickerson-Drew DNA is comparable to the scatter of 3JP,H3' couplings among four separate NMR experiments. A commonly accepted assumption of homogeneity of one kind of structure-dynamic state within DNA (BI or BII) is questionable because the principal characteristics of relevant P(ε) probabilities (shapes and positioning) vary with DNA sequence. The theory outlined in the present work sets limits to future reparameterization of MD force fields, as relevant to NMR data.

The Ad-MD method to calculate NMR shift including effects due to conformational dynamics: The 31 P NMR shift in DNA.

A method for averaging of NMR parameters by molecular dynamics (MD) has been derived from the method of statistical averaging in MD snapshots, benchmarked and applied to structurally dynamic interpretation of the 31 P NMR shift (δ31P ) in DNA phosphates. The method employs adiabatic dependence of an NMR parameter on selected geometric parameter(s) that is weighted by MD-calculated probability distribution(s) for the geometric parameter(s) (Ad-MD method). The usage of Ad-MD for polymers is computationally convenient when one pre-calculated structural dependence of an NMR parameter is employed for all chemically equivalent units differing only in dynamic behavior. The Ad-MD method is benchmarked against the statistical averaging method for δ31P in the model phosphates featuring distinctively different structures and dynamic behavior. The applicability of Ad-MD is illustrated by calculating 31 P NMR spectra in the Dickerson-Drew DNA dodecamer. δ31P was calculated with the B3LYP/IGLO-III/PCM(water) and the probability distributions for the torsion angles adjacent to the phosphorus atoms in the DNA phosphates were calculated using the OL15 force field.

Z-DNA as a Touchstone for Additive Empirical Force Fields and a Refinement of the Alpha/Gamma DNA Torsions for AMBER.

Although current AMBER force fields are relatively accurate for canonical B-DNA, many noncanonical structures are still described incorrectly. As noncanonical motifs are attracting increasing attention due to the role they play in living organisms, further improvement is desirable. Here, we have chosen the Z-DNA molecule, which can be considered a touchstone of the universality of empirical force fields, since the noncanonical α and γ backbone conformations native to Z-DNA are also found in protein-DNA complexes, i-motif DNA, and other noncanonical DNAs. We show that spurious α/γ conformations occurring in simulations with current AMBER force fields, OL15 and bsc1, are largely due to inaccurate α/γ parametrization. Moreover, stabilization of native Z-DNA substates involving γ = trans conformations appears to be in conflict with the correct description of the canonical B-DNA structure. Because the balance of the native and spurious conformations is influenced by nonadditive effects, this is a difficult case for an additive dihedral energy scheme such as AMBER. We propose new α/γ parameters, denoted OL21, and show that they improve the stability of native α/γ Z-DNA substates while keeping the canonical DNA description virtually unchanged, thus representing a reasonable compromise within the additive force field framework. Although further extensive testing is needed, the new modification appears to be a promising step toward a more reliable description of noncanonical DNA motifs and provides the best performance for Z-DNA molecules among current AMBER force fields.

Improving the Performance of the Amber RNA Force Field by Tuning the Hydrogen-Bonding Interactions.

Molecular dynamics (MD) simulations became a leading tool for investigation of structural dynamics of nucleic acids. Despite recent efforts to improve the empirical potentials (force fields, ffs), RNA ffs have persisting deficiencies, which hamper their utilization in quantitatively accurate simulations. Previous studies have shown that at least two salient problems contribute to difficulties in the description of free-energy landscapes of small RNA motifs: (i) excessive stabilization of the unfolded single-stranded RNA ensemble by intramolecular base-phosphate and sugar-phosphate interactions and (ii) destabilization of the native folded state by underestimation of stability of base pairing. Here, we introduce a general ff term (gHBfix) that can selectively fine-tune nonbonding interaction terms in RNA ffs, in particular, the H bonds. The gHBfix potential affects the pairwise interactions between all possible pairs of the specific atom types, while all other interactions remain intact; i.e., it is not a structure-based model. In order to probe the ability of the gHBfix potential to refine the ff nonbonded terms, we performed an extensive set of folding simulations of RNA tetranucleotides and tetraloops. On the basis of these data, we propose particular gHBfix parameters to modify the AMBER RNA ff. The suggested parametrization significantly improves the agreement between experimental data and the simulation conformational ensembles, although our current ff version still remains far from being flawless. While attempts to tune the RNA ffs by conventional reparametrizations of dihedral potentials or nonbonded terms can lead to major undesired side effects, as we demonstrate for some recently published ffs, gHBfix has a clear promising potential to improve the ff performance while avoiding introduction of major new imbalances.

A- to B-DNA Transition in AMBER Force Fields and Its Coupling to Sugar Pucker.

The A/B transition is a basic element of DNA conformational change. Because of its involvement in the sensing of the ionic conditions by DNA and in specific protein-DNA interactions, this transition is important for biological functions of DNA. Therefore, accurate modeling of the A/B equilibrium by means of empirical force fields is of utmost interest. In this work, we examine the A/B equilibrium in three AMBER force fields, including the recent bsc1 and OL15 modifications, using much longer MD simulations than attempted before. Special attention is paid to the coupling of the A/B equilibrium with the south/north (S/N) transition of the sugar pucker. We found that none of the tested force fields provided a satisfactory description of the A/B equilibrium because the B-form was predicted to be much too stable and the A-form was predicted to be almost absent even in concentrated trifluoroethanol solutions. Based on comparison with NMR data for duplexes and single nucleosides, we hypothesize that this problem arose from the incorrect description of the S/N equilibrium of sugar pucker, where the south conformation is much too stable, thus stabilizing the B-form. Because neither the A/B equilibrium in duplexes nor the S/N equilibrium in nucleosides was described accurately, further refinements of the AMBER DNA force fields are needed.

Noncanonical α/γ Backbone Conformations in RNA and the Accuracy of Their Description by the AMBER Force Field.

The sugar-phosphate backbone of RNA can exist in diverse rotameric substates, giving RNA molecules enormous conformational variability. The most frequent noncanonical backbone conformation in RNA is α/γ = t/t, which is derived from the canonical backbone by a crankshaft motion and largely preserves the standard geometry of the RNA duplex. A similar conformation also exists in DNA, where it has been extensively studied and shown to be involved in DNA-protein interactions. However, the function of the α/γ = t/t conformation in RNA is poorly understood. Here, we present molecular dynamics simulations of several prototypical RNA structures obtained from X-ray and NMR experiments, including canonical and mismatched RNA duplexes, UUCG and GAGA tetraloops, Loop E, the sarcin-ricin loop, a parallel guanine quadruplex, and a viral pseudoknot. The stability of various noncanonical α/γ backbone conformations was analyzed with two AMBER force fields, ff99bsc0χOL3 and ff99bsc0χOL3 with the recent εζOL1 and ßOL1 corrections for DNA. Although some α/γ substates were stable with seemingly well-described equilibria, many were unstable in our simulations. Notably, the most frequent noncanonical conformer α/γ = t/t was unstable in both tested force fields. Possible reasons for this instability are discussed. Our work reveals a potentially important artifact in RNA force fields and highlights a need for further force field refinement.

Influence of BII Backbone Substates on DNA Twist: A Unified View and Comparison of Simulation and Experiment for All 136 Distinct Tetranucleotide Sequences.

Reliable representation of the B-DNA base-pair step twist is one of the crucial requirements for theoretical modeling of DNA supercoiling and other biologically relevant phenomena in B-DNA. It has long been suspected that the twist is inaccurately described by current empirical force fields. Unfortunately, comparison of simulation results with experiments is not straightforward because of the presence of BII backbone substates, whose populations may differ in experimental and simulation ensembles. In this work, we provide a comprehensive view of the effect of BII substates on the overall B-DNA helix twist and show how to reliably compare twist values from experiment and simulation in two scenarios. First, for longer DNA segments freely moving in solution, we show that sequence-averaged twists of different BI/BII ensembles can be compared directly because of approximate cancellation of the opposing BII effects. Second, for sequence-specific data, such as a particular base-pair step or tetranucleotide twist, can be compared only for a clearly defined BI/BII backbone conformation. For the purpose of force field testing, we designed a compact set of fourteen 22-base-pair B-DNA duplexes (Set 14) containing all 136 distinct tetranucleotide sequences and carried out a total of 84 µs of molecular dynamics simulations, primarily with the OL15 force field. Our results show that the ff99bsc0εζOL1χOL4, parmbsc1, and OL15 force fields model the B-DNA helical twist in good agreement with X-ray and minicircle ligation experiments. The comprehensive understanding obtained regarding the effect of BII substates on the base-pair step geometry should aid meaningful comparisons of various conformational ensembles in future research.

How to understand atomistic molecular dynamics simulations of RNA and protein-RNA complexes?

We provide a critical assessment of explicit-solvent atomistic molecular dynamics (MD) simulations of RNA and protein/RNA complexes, written primarily for non-specialists with an emphasis to explain the limitations of MD. MD simulations can be likened to hypothetical single-molecule experiments starting from single atomistic conformations and investigating genuine thermal sampling of the biomolecules. The main advantage of MD is the unlimited temporal and spatial resolution of positions of all atoms in the simulated systems. Fundamental limitations are the short physical time-scale of simulations, which can be partially alleviated by enhanced-sampling techniques, and the highly approximate atomistic force fields describing the simulated molecules. The applicability and present limitations of MD are demonstrated on studies of tetranucleotides, tetraloops, ribozymes, riboswitches and protein/RNA complexes. Wisely applied simulations respecting the approximations of the model can successfully complement structural and biochemical experiments. WIREs RNA 2017, 8:e1405. doi: 10.1002/wrna.1405 For further resources related to this article, please visit the WIREs website.

Assessing the Current State of Amber Force Field Modifications for DNA.

The utility of molecular dynamics (MD) simulations to model biomolecular structure, dynamics, and interactions has witnessed enormous advances in recent years due to the availability of optimized MD software and access to significant computational power, including GPU multicore computing engines and other specialized hardware. This has led researchers to routinely extend conformational sampling times to the microsecond level and beyond. The extended sampling time has allowed the community not only to converge conformational ensembles through complete sampling but also to discover deficiencies and overcome problems with the force fields. Accuracy of the force fields is a key component, along with sampling, toward being able to generate accurate and stable structures of biopolymers. The Amber force field for nucleic acids has been used extensively since the 1990s, and multiple artifacts have been discovered, corrected, and reassessed by different research groups. We present a direct comparison of two of the most recent and state-of-the-art Amber force field modifications, bsc1 and OL15, that focus on accurate modeling of double-stranded DNA. After extensive MD simulations with five test cases and two different water models, we conclude that both modifications are a remarkable improvement over the previous bsc0 force field. Both force field modifications show better agreement when compared to experimental structures. To ensure convergence, the Drew-Dickerson dodecamer (DDD) system was simulated using 100 independent MD simulations, each extended to at least 10 µs, and the independent MD simulations were concatenated into a single 1 ms long trajectory for each combination of force field and water model. This is significantly beyond the time scale needed to converge the conformational ensemble of the internal portions of a DNA helix absent internal base pair opening. Considering all of the simulations discussed in the current work, the MD simulations performed to assess and validate the current force fields and water models aggregate over 14 ms of simulation time. The results suggest that both the bsc1 and OL15 force fields render average structures that deviate significantly less than 1 Å from the average experimental structures. This can be compared to similar but less exhaustive simulations with the CHARMM 36 force field that aggregate to the ∼90 µs time scale and also perform well but do not produce structures as close to the DDD NMR average structures (with root-mean-square deviations of 1.3 Å) as the newer Amber force fields. On the basis of these analyses, any future research involving double-stranded DNA simulations using the Amber force fields should employ the bsc1 or OL15 modification.

On the Use of Molecular Dynamics Simulations for Probing Allostery through DNA.

A recent study described an allosteric effect in which the binding of a protein to DNA is influenced by another protein bound nearby. The effect shows a periodicity of ∼10 basepairs and decays with increasing protein-protein distance. As a mechanistic explanation, the authors reported a similar periodic, decaying pattern of the correlation coefficient between major groove widths inferred from a shorter molecular dynamics simulation. Here we show that in a state-of-the-art, microsecond-long simulation of the same DNA sequence, the periodicity of the correlation coefficient is not observed. To study the problem further, we extend an earlier mechanical model of DNA allostery based on constrained minimization of effective quadratic deformation energy of the DNA. We demonstrate that, if the constraints mimicking the bound proteins are properly applied, the periodicity in the binding energy is indeed recovered.

Refinement of the Sugar-Phosphate Backbone Torsion Beta for AMBER Force Fields Improves the Description of Z- and B-DNA.

Z-DNA duplexes are a particularly complicated test case for current force fields. We performed a set of explicit solvent molecular dynamics (MD) simulations with various AMBER force field parametrizations including our recent refinements of the ε/ζ and glycosidic torsions. None of these force fields described the ZI/ZII and other backbone substates correctly, and all of them underpredicted the population of the important ZI substate. We show that this underprediction can be attributed to an inaccurate potential for the sugar-phosphate backbone torsion angle ß. We suggest a refinement of this potential, ß(OL1), which was derived using our recently introduced methodology that includes conformation-dependent solvation effects. The new potential significantly increases the stability of the dominant ZI backbone substate and improves the overall description of the Z-DNA backbone. It also has a positive (albeit small) impact on another important DNA form, the antiparallel guanine quadruplex (G-DNA), and improves the description of the canonical B-DNA backbone by increasing the population of BII backbone substates, providing a better agreement with experiment. We recommend using ß(OL1) in combination with our previously introduced corrections, εζ(OL1) and χ(OL4), (the combination being named OL15) as a possible alternative to the current ß torsion potential for more accurate modeling of nucleic acids.

Microsecond-Scale MD Simulations of HIV-1 DIS Kissing-Loop Complexes Predict Bulged-In Conformation of the Bulged Bases and Reveal Interesting Differences between Available Variants of the AMBER RNA Force Fields.

We report an extensive set of explicit solvent molecular dynamics (MD) simulations (∼25 µs of accumulated simulation time) of the RNA kissing-loop complex of the HIV-1 virus initiation dimerization site. Despite many structural investigations by X-ray, NMR, and MD techniques, the position of the bulged purines of the kissing complex has not been unambiguously resolved. The X-ray structures consistently show bulged-out positions of the unpaired bases, while several NMR studies show bulged-in conformations. The NMR studies are, however, mutually inconsistent regarding the exact orientations of the bases. The earlier simulation studies predicted the bulged-out conformation; however, this finding could have been biased by the short simulation time scales. Our microsecond-long simulations reveal that all unpaired bases of the kissing-loop complex stay preferably in the interior of the kissing-loop complex. The MD results are discussed in the context of the available experimental data and we suggest that both conformations are biochemically relevant. We also show that MD provides a quite satisfactory description of this RNA system, contrasting recent reports of unsatisfactory performance of the RNA force fields for smaller systems such as tetranucleotides and tetraloops. We explain this by the fact that the kissing complex is primarily stabilized by an extensive network of Watson-Crick interactions which are rather well described by the force fields. We tested several different sets of water/ion parameters but they all lead to consistent results. However, we demonstrate that a recently suggested modification of van der Waals interactions of the Cornell et al. force field deteriorates the description of the kissing complex by the loss of key stacking interactions stabilizing the interhelical junction and excessive hydrogen-bonding interactions.

Transferability and additivity of dihedral parameters in polarizable and nonpolarizable empirical force fields.

Recent advances in polarizable force fields have revealed that major reparameterization is necessary when the polarization energy is treated explicitly. This study is focused on the torsional parameters, which are crucial for the accurate description of conformational equilibria in biomolecules. In particular, attention is paid to the influence of polarization on the (i) transferability of dihedral terms between molecules, (ii) transferability between different environments, and (iii) additivity of dihedral energies. To this end, three polarizable force fields based on the induced point dipole model designed for use in AMBER are tested, including two recent ff02 reparameterizations. Attention is paid to the contributions due to short range interactions (1-2, 1-3, and 1-4) within the four atoms defining the dihedral angle. The results show that when short range 1-2 and 1-3 polarization interactions are omitted, as for instance in ff02, the 1-4 polarization contribution is rather small and unlikely to improve the description of the torsional energy. Conversely, when screened 1-2 and 1-3 interactions are included, the polarization contribution is sizeable and shows potential to improve the transferability of parameters between different molecules and environments as well as the additivity of dihedral terms. However, to reproduce intramolecular polarization effects accurately, further fine-tuning of the short range damping of polarization is necessary.

Reactive conformation of the active site in the hairpin ribozyme achieved by molecular dynamics simulations with ε/ζ force field reparametrizations.

X-ray crystallography can provide important insights into the structure of RNA enzymes (ribozymes). However, the details of a ribozyme's active site architecture are often altered by the inactivating chemical modifications necessary to inhibit self-cleavage. Molecular dynamics (MD) simulations are able to complement crystallographic data and model the conformation of the ribozyme's active site in its native form. However, the performance of MD simulations is driven by the quality of the force field used. Force fields are primarily parametrized and tested for a description of canonical structures and thus may be less accurate for noncanonical RNA elements, including ribozyme catalytic cores. Here, we show that our recent reparametrization of ε/ζ torsions significantly improves the description of the hairpin ribozyme's scissile phosphate conformational behavior. In addition, we find that an imbalance in the force field description of the nonbonded interactions of the ribose 2'-OH contributes to the conformational behavior observed for the scissile phosphate in the presence of a deprotonated G8(-). On the basis of the new force field, we obtain a reactive conformation for the hairpin ribozyme active site that is consistent with the most recent mechanistic and structural data.

Mechanical properties of symmetric and asymmetric DNA A-tracts: implications for looping and nucleosome positioning.

A-tracts are functionally important DNA sequences which induce helix bending and have peculiar structural properties. While A-tract structure has been qualitatively well characterized, their mechanical properties remain controversial. A-tracts appear structurally rigid and resist nucleosome formation, but seem flexible in DNA looping. In this work, we investigate mechanical properties of symmetric AnTn and asymmetric A2n tracts for n = 3, 4, 5 using two types of coarse-grained models. The first model represents DNA as an ensemble of interacting rigid bases with non-local quadratic deformation energy, the second one treats DNA as an anisotropically bendable and twistable elastic rod. Parameters for both models are inferred from microsecond long, atomic-resolution molecular dynamics simulations. We find that asymmetric A-tracts are more rigid than the control G/C-rich sequence in localized distortions relevant for nucleosome formation, but are more flexible in global bending and twisting relevant for looping. The symmetric tracts, in contrast, are more rigid than asymmetric tracts and the control, both locally and globally. Our results can reconcile the contradictory stiffness data on A-tracts and suggest symmetric A-tracts to be more efficient in nucleosome exclusion than the asymmetric ones. This would open a new possibility of gene expression manipulation using A-tracts.

Molecular Dynamics Simulations of Nucleic Acids. From Tetranucleotides to the Ribosome.

We present a brief overview of explicit solvent molecular dynamics (MD) simulations of nucleic acids. We explain physical chemistry limitations of the simulations, namely, the molecular mechanics (MM) force field (FF) approximation and limited time scale. Further, we discuss relations and differences between simulations and experiments, compare standard and enhanced sampling simulations, discuss the role of starting structures, comment on different versions of nucleic acid FFs, and relate MM computations with contemporary quantum chemistry. Despite its limitations, we show that MD is a powerful technique for studying the structural dynamics of nucleic acids with a fast growing potential that substantially complements experimental results and aids their interpretation.

Mechanical Model of DNA Allostery.

The importance of allosteric effects in DNA is becoming increasingly appreciated, but the underlying mechanisms remain poorly understood. In this work, we propose a general modeling framework to study DNA allostery. We describe DNA in a coarse-grained manner by intra-base pair and base pair step coordinates, complemented by groove widths. Quadratic deformation energy is assumed, yielding linear relations between the constraints and their effect. Model parameters are inferred from standard unrestrained, explicit-solvent molecular dynamics simulations of naked DNA. We applied the approach to study minor groove binding of diamidines and pyrrole-imidazole polyamides. The predicted DNA bending is in quantitative agreement with experiment and suggests that diamidine binding to the alternating TA sequence brings the DNA closer to the A-tract conformation, with potentially important functional consequences. The approach can be readily applied to other allosteric effects in DNA and generalized to model allostery in various molecular systems.